% This code produces the main figure for the monetary part of the paper
% but tries out different monetary policy rules
% Here: we increase phi_pi in the Taylor rule to large numbers
% by Raphael schoenle, February 1, 2016
% edited by Jiahceng Feng modified March 18, 2019

% six cases with parameterizations

clearvars -except wd

% Data
load FPA.mat
%Parameters values
params.rho_a = 0.9;         % persistence agg tech shocks
params.rho_mu = 0.9;         % persistence monetary shock
params.rho_k = 0.9;         % persistence idiosync tech shock
params.beta = 0.9975;      % discount factor
params.theta = 6;          % elast of subst within sectors
params.delta = 0.5;

Ck = csvread("2002CkDetailed.csv");
Znum = csvread("2002RevShareDetailed.csv");
n = length(Ck); % Number of Industries
Z = Znum(1:n,1:n);

Wc = Ck/sum(Ck); % Share of consumption; the weights omega
W = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega
psi = params.delta*(params.theta - 1)/params.theta;
psi_C = (1-params.delta*(params.theta - 1)/params.theta);
n_k = inv(eye(n)-psi.*W)*(psi_C*Wc);

%number of sectors
n_sectors = 341;
n_sectors1 = 56;

%number of periods in IRF
n_periods = 30;
n_periods_idio = 30;

%number of cases
num_case = 18;

%pre-definition of variables
% third dimension: number of shocks: monetary, technology, idiosyncratic
c_imp = zeros(n_periods, num_case, 3);
inf_imp = zeros(n_periods, num_case, 3);
mc_imp = zeros(n_periods, num_case, 3, n_sectors);
c_imp_old=c_imp;
inf_imp_old=inf_imp;
mc_imp_old=mc_imp;
pos_akshock = zeros(n_sectors,1);
pos_akshock1 = zeros(n_sectors1,1);

% parameter sets
delta_set = [ 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
eta_set = [ 2 2 2 2 2 2 2];
phi_set = [2 2 2 2 2 2 2];
sigma_set = [ 1 1 1 1 1 1 1];
phi_pi_set = [1.24 1.24 1.24 1.24 1.24 1.24 1.24];
phi_c_set = [ 0.33 0.33 0.33 0.33 0.33 0.33 0.33]./12;
phi_rho_i_set = [ 0 0 0.1 0.5 0.7 0.8 0.9];
phi_gc_i_set = [ 0 1.5 1.5 1.5 1.5 1.5 1.5];
n_parameters = length(delta_set);

% Position indices for shocks
pos_mushock = 1;
pos_ashock = 2;
for i = 1:n_sectors
    pos_akshock(i) = 2 + i;
end
for i = 1:n_sectors1
    pos_akshock1(i) = 2 + i;
end

% sectoral weights in consumption
Omega_C_het = Wc;
Omega_C_hom = 1/n_sectors*ones(size(Omega_C_het));

% load sector code labels
% initial data
 agg_1 = fopen('data_341_NAICS_sector_names.csv');
 fmt = '%s%s';
 data_1 = textscan(agg_1,fmt,'headerlines',1,'delimiter',',');
% BEA NAICS codes: these are always 341 identical codes from the first column
 bea_codes_341 = data_1{1};

% sectoral Calvo
% 341 sectors first
FPA_het  = FPA;
FPA_hom = mean(FPA_het)*ones(size(FPA_het)); 
FPA_flex = 0.99*ones(size(FPA_het));
FPA_sticky = 0.01*ones(size(FPA_het));
% Omega_C weighted mean
FPA_hom_weighted = Omega_C_het'*FPA_het*ones(size(FPA_het));

%sectoral weights 
Omega_het = W;
Omega_hom = 1/n_sectors*ones(size(Omega_het)); 
Omega_C_as_IO = repmat(Omega_C_het',n_sectors,1);

%% 341 Sectors
%Case1: flex-price Calvo, homogeneous size, homogeneous I/O
FPA(:,:,1) = FPA_flex;
Omega_C(:,:,1) = Omega_C_hom;
Omega(:,:,1) = Omega_hom;

%Case2: homogeneous weighted mean Calvo, homogeneous size, homogeneous I/O
FPA(:,:,2) = FPA_hom_weighted;
Omega_C(:,:,2) = Omega_C_hom;
Omega(:,:,2) = Omega_hom;

%Case3: heterogeneous Calvo, homogeneous size, homogeneous I/O
FPA(:,:,3) = FPA_het;
Omega_C(:,:,3) = Omega_C_hom;
Omega(:,:,3) = Omega_hom;

%Case4: heterogeneous Calvo, heterogeneous size, heterogeneous
%I/O based on size weights
FPA(:,:,4) = FPA_het;
Omega_C(:,:,4) = Omega_C_het;
Omega(:,:,4) = Omega_C_as_IO;

%Case5: heterogeneous Calvo, heterogeneous size, homogeneous I/O
FPA(:,:,5) = FPA_het;
Omega_C(:,:,5) = Omega_C_het;
Omega(:,:,5) = Omega_hom;

%Case6: heterogeneous Calvo, heterogeneous size, heterogeneous I/O
FPA(:,:,6) = FPA_het;
Omega_C(:,:,6) = Omega_C_het;
Omega(:,:,6) = Omega_het;

% outer loop over the parameter values, index k
for k = 6:n_parameters
% baseline
%for k = 2:2
    
    params.delta = delta_set(k);
    params.eta = eta_set(k);
    params.frisch = phi_set(k);
    params.sigma = sigma_set(k);
    params.phi_pi = phi_pi_set(k);
    params.phi_c = phi_c_set(k);
    params.phi_gc = phi_gc_i_set(k);
    params.rho_i = phi_rho_i_set(k);

for i = 6:6

        % delete Table 5a,b output if existing
        if exist(['Table5a_case_' num2str(i) '_params' num2str(k) '_341.csv'], 'file')==2
            delete(['Table5a_case_' num2str(i) '_params' num2str(k) '_341.csv']);
        end
        if exist(['Table5b_case_' num2str(i) '_params' num2str(k) '_341.csv'], 'file')==2
            delete(['Table5b_case_' num2str(i) '_params' num2str(k) '_341.csv']);
        end
        tic

        % Generate Gensys matrices
        [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps]  = ...
                       Gensys_matrix(n_sectors, FPA(:,:,i),Omega_C(:,:,i), Omega(:,:,i), params);
%            Gensys_matrix_interest_rate_smoothing(n_sectors, FPA(:,:,i),Omega_C(:,:,i), Omega(:,:,i), params);
                      
    toc
            %%Call Gensys
            [G1, C, imp, fmat, fwt, ywt, gev, eu, loose] = gensys(G0, G1, C, PSI, PI);
        toc
            %check existence and uniqueness
            eu
    
        %%Generate IRFs
        %%%Monetary policy shocks
        % now include MC as output
            yt_imp_mu = Gensys_IRF(G1, n_sectors, n_periods, pos_mushock, imp);
            % output = [index for sectors, consumption IRF, inflation IRF, real MC IRF]
            output = [1:n_periods; yt_imp_mu(pos_c,2:end); ...
                yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1); ...
                yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)) ]';
            autocorr_c = autocorr(output(:,end-1),1);
            autocorr_inf = autocorr(output(:,end),1);
            autocorr_mc = autocorr([mean(yt_imp_mu(pos_mc,2:end)) - (yt_imp_mu(pos_pc,2:end))]',1);

            % output matrix: initial C, Pi, MC, cumulative C, Pi, mean(MC)...
            % AR(1) of C, Pi, MC + IRFs (C, Pi, MC)
            output_matrix = [output(:,1) repmat([output(1,2:3)...
                mean(yt_imp_mu(pos_mc,2) - (yt_imp_mu(pos_pc,2)))...
            sum(output(:,2:3)) sum(yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)))...    
            autocorr_c(2) autocorr_inf(2) autocorr_mc(2)],n_periods,1) output(:,2:end)];

            % cumulative responses - sort and get sector indices
            [output_top_bottom, I] = sort(sum(yt_imp_mu(pos_yk,2:end),2));
            top_bottom = [I(1:10) ; I(end-9:end)];

        % save cumulative response and identity of extreme bottom/top 10 sectors
        % output all sectors 
        filename = ['Table5a_case_' num2str(i) 'params' num2str(k) '_341.csv'];
        dlmwrite(filename, [output_top_bottom]', '-append',  'precision', '%.10f', 'delimiter', '\t');
%        dlmwrite(filename, [output_top_bottom(1:10) output_top_bottom(end-9:end)]', '-append',  'precision', '%.10f', 'delimiter', '\t');
        [output_top_bottom(1:10) output_top_bottom(end-9:end)]'

        filename = ['Table5b_case_' num2str(i) 'params' num2str(k) '_341.csv'];
        fid = fopen(filename, 'a');
        %fprintf(fid, '"%s"\n', [ strjoin(bea_codes_341(I))]);
        fclose(fid)
        
        % save output
        filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.csv'];
        fid = fopen(filename, 'w');
        fprintf(fid,'Period\tConsumption_impact_irf_mu\tInflation_impact_irf_mu\tMC_mean_impact_irf_mu\tConsumption_cum_irf_mu\tInflation_cum_irf_mu\tMC_mean_cum_irf_mu\tConsumption_persistence_mu\tInflation_persistence_mu\tMC_persistence_mu\tConsumption_mu\tInflation_mu\tMC_mu \n');
        fclose(fid);
        dlmwrite(filename, (output_matrix), '-append', 'precision', '%.10f', 'delimiter', '\t');
        
        filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat']
        save(filename, 'output_matrix')
       
        % for plotting
        c_imp(:,i,1) = yt_imp_mu(pos_c,2:end)';
        inf_imp(:,i,1) = (yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1))';
        mc_imp(:,i,1) = [mean(yt_imp_mu(pos_mc,2:end)) - yt_imp_mu(pos_pc,2:end)]';
 
    toc
    end
    
end

 
save('alloutput_main_interest_rate_smoothing.mat')

% we plot baseline, and and acses 3, 5, 7

%cd '/Users/raphaelschoenle/Dropbox/Networks/code/JME_post/Figure2'

load alloutput_main_interest_rate_smoothing.mat

% parameter sets
delta_set = [ 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
eta_set = [ 2 2 2 2 2 2 2];
phi_set = [2 2 2 2 2 2 2];
sigma_set = [ 1 1 1 1 1 1 1];
phi_pi_set = [1.24 1.24 1.24 1.24 1.24 1.24 1.24];
phi_c_set = [ 0.33 0.33 0.33 0.33 0.33 0.33 0.33]./12;
phi_rho_i_set = [ 0 0 0.1 0.5 0.7 0.8 0.9];
phi_gc_i_set = [ 0 1.5 1.5 1.5 1.5 1.5 1.5];

i=6

% case 1/ baseline 
k=1
filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
load(filename)    
n_sectors = 341
c1 = output_matrix(:,end-n_sectors-1);
pi1 = output_matrix(:,end-n_sectors);
mc1 =  mean(output_matrix(:,end-n_sectors+1:end),2);

% interest rate smoothing 0.1
k=3
filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
load(filename)   
c2 = output_matrix(:,end-n_sectors-1);
pi2 = output_matrix(:,end-n_sectors);
mc2 =  mean(output_matrix(:,end-n_sectors+1:end),2);

% interest rate smoothing 0.7
k=5
filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
load(filename)   
c3 = output_matrix(:,end-n_sectors-1);
pi3 = output_matrix(:,end-n_sectors);
mc3 =  mean(output_matrix(:,end-n_sectors+1:end),2);

% interest rate smoothing 0.9
k=7
filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
load(filename)   
c4 = output_matrix(:,end-n_sectors-1);
pi4 = output_matrix(:,end-n_sectors);
mc4 = mean(output_matrix(:,end-n_sectors+1:end),2);
  
% % Parameters used for MP rule:
% %phi_pi_set = [1.24 1.24 1.24 1.24 1.24 1.24];
% %phi_c_set = [0 0.33 0.33 0.33 0.33 0.33]./12;
% %phi_rho_i_set = [0 0 0 0.1 0.5 0.9];
% %phi_gc_i_set = [0 0 1.5 1.5 1.5 1.5];
% 

% we adjust the initial impact responses for interest rate smoothing by
% choosing a shock size such that we get the same initial response as for
% flexible inflation targeting (but not including output growth)
figure('units','normalized','position',[.1 0 .4 1])
subplot(3,1,1)
plot(c1,'-b','linewidth',2)
hold on
plot(c2*c1(1,1)/c2(1,1),'-b+','linewidth',2)
plot(c3*c1(1,1)/c3(1,1),'--b','linewidth',2)
plot(c4*c1(1,1)/c4(1,1),':b','linewidth',2)
    title('\textbf{Consumption Response}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')

subplot(3,1,2)
plot(pi1,'-b','linewidth',2)
hold on
plot(pi2*c1(1,1)/c2(1,1),'-b+','linewidth',2)
plot(pi3*c1(1,1)/c3(1,1),'--b','linewidth',2)
plot(pi4*c1(1,1)/c4(1,1),':b','linewidth',2)
    legend1=legend('\phi_{\pi}=1.24 \phi_{c}=0.33/12',...
        '\phi_{\pi}=1.24 \phi_{c}=0.33/12 \phi_{gc}=1.5 \rho_i=0.1',...
        '\phi_{\pi}=1.24 \phi_{c}=0.33/12 \phi_{gc}=1.5 \rho_i=0.5',...
        '\phi_{\pi}=1.24 \phi_{c}=0.33/12 \phi_{gc}=1.5 \rho_i=0.9',.....
        'Location','best')
    title('\textbf{Inflation Response}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')

subplot(3,1,3)
plot(mc1,'-b','linewidth',2)
hold on
plot(mc2*c1(1,1)/c2(1,1),'-b+','linewidth',2)
plot(mc3*c1(1,1)/c3(1,1),'--b','linewidth',2)
plot(mc4*c1(1,1)/c4(1,1),':b','linewidth',2)
    title('\textbf{Real MC}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    fig341 = gcf
    fig341.PaperUnits = 'normalized';
    fig341.PaperPosition = [.1 0 .8 1];
    fig341.PaperPositionMode = 'manual';
     filename = ['IRFS_full_heterogeneity_MP_rules_normalized_smoothing', '.pdf'];
     saveas(fig341,filename);
     filename = ['IRFS_full_heterogeneity_MP_rules_normalized_smoothing', '.eps'];
     saveas(fig341,filename);
